Deep dive regional communities - using relative abundance prox

source('../env.R')

Species in communities

It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.

Load data

city_effort = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'city_effort.csv'))

Rows: 342 Columns: 7── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): city_name
dbl (6): city_id, total_city_checklists, total_city_locations, total_city_effort, total_city_area_m2, percentage_total_city_area_surveyed
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

city_effort

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2462 Columns: 12── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (4): city_id, distance_to_northern_edge_km, distance_to_southern_edge_km, relative_abundance_proxy
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities_summary = communities %>% group_by(city_id) %>% summarise(
  regional_pool_size = n(), 
  urban_pool_size = sum(relative_abundance_proxy > 0)
) %>% left_join(city_effort %>% dplyr::select(city_id, percentage_total_city_area_surveyed))

Joining with `by = join_by(city_id)`

ggplot(communities %>% filter(relative_abundance_proxy > 0), aes(x = relative_abundance_proxy)) + geom_bar(stat = "bin")

city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))

Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude data

community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv')) %>%
  dplyr::select(city_id, mntd_normalised, fdiv_normalised, mass_fdiv_normalised, handwing_index_fdiv_normalised, trophic_trait_fdiv_normalised, gape_width_fdiv_normalised) %>%
  left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
  left_join(communities_summary) %>%
  left_join(city_points[,c('city_id', 'city_nm')]) %>%
  rename(city_name='city_nm') %>%
  na.omit() %>%
  arrange(city_id)

Rows: 341 Columns: 43── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (43): mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fdiv_normalised, fdiv_actual, fdiv_min, fdiv_max, fdiv_mean, fdiv_sd, mass_fdi...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`

community_data_metrics

community_data_metrics %>% group_by(core_realm) %>% summarise(total_cities = n())

test_value = function(name, normalised_list) {
  wilcox_test_result = wilcox.test(normalised_list, mu = 0.5)
  
  significance = ifelse(wilcox_test_result$p.value < 0.0001, '***', 
                        ifelse(wilcox_test_result$p.value < 0.001, '**', 
                               ifelse(wilcox_test_result$p.value < 0.01, '*', '')))
  m = mean(normalised_list)
  
  paste(name, 'mean', round(m, 2), significance)
}

test_value('Global MNTD', community_data_metrics$mntd_normalised)

[1] "Global MNTD mean 0.5 "

test_value('Global beak gape FDiv', community_data_metrics$gape_width_fdiv_normalised)

[1] "Global beak gape FDiv mean 0.59 ***"

test_value('Global HWI FDiv', community_data_metrics$handwing_index_fdiv_normalised)

[1] "Global HWI FDiv mean 0.61 ***"

nrow(community_data_metrics)

[1] 319

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))

Rows: 304 Columns: 6── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (5): gape_width, trophic_trait, locomotory_trait, mass, handwing_index
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

fetch_normalised_traits = function(required_species_list) {
  required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
  
  required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
  required_traits$handwing_index_fdiv_normalised = normalise(required_traits$handwing_index, min(required_traits$handwing_index), max(required_traits$handwing_index))
  required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
  
  traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'handwing_index_fdiv_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
  traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'handwing_index_fdiv_normalised', 'mass_normalised'), labels = c('Gape Width', 'HWI', 'Mass'))
  
  traits_normalised_long
}

fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))

Read in our phylogeny

phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')

Load resolve ecoregions

resolve = read_resolve()

Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

Create helper functions

to_species_matrix = function(filtered_communities) {
  filtered_communities %>% 
    dplyr::select(city_id, jetz_species_name) %>% 
    distinct() %>%
    mutate(present = TRUE) %>% 
    pivot_wider(
      names_from = jetz_species_name, 
      values_from = "present", 
      values_fill = list(present = F)
    ) %>% 
    tibble::column_to_rownames(var='city_id')
}

community_nmds = function(filtered_communities) {
  species_matrix = to_species_matrix(filtered_communities)
  nmds <- metaMDS(species_matrix, k=2, trymax = 30)
  nmds_result = data.frame(vegan::scores(nmds)$sites)
  nmds_result$city_id = as.double(rownames(nmds_result))
  rownames(nmds_result) = NULL
  nmds_result
}

https://www.datacamp.com/tutorial/k-means-clustering-r

scree_plot = function(community_nmds_data) {
  # Decide how many clusters to look at
  n_clusters <- min(10, nrow(community_nmds_data) - 1)
  
  # Initialize total within sum of squares error: wss
  wss <- numeric(n_clusters)
  
  set.seed(123)
  
  # Look over 1 to n possible clusters
  for (i in 1:n_clusters) {
    # Fit the model: km.out
    km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
    # Save the within cluster sum of squares
    wss[i] <- km.out$tot.withinss
  }
  
  # Produce a scree plot
  wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
   
  scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
      geom_point(size = 4) +
      geom_line() +
      geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
      scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
      xlab('Number of clusters')
  scree_plot
}

cluster_cities = function(city_nmds, cities_community_data, centers) {
  set.seed(123)
  kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
  city_nmds$cluster = kmeans_clusters$cluster
  cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}

plot_nmds_clusters = function(cluster_cities) {
  cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}

get_presence_cell_width = function(city_cluster_data_metrics) {
  10 * length(unique(city_cluster_data_metrics$city_id))
}

get_presence_cell_height = function(city_cluster_data_metrics) {
  species = species_in_cluster = communities %>% 
    filter(city_id %in% city_cluster_data_metrics$city_id) %>% 
    dplyr::select(jetz_species_name) %>% 
    distinct()
  
  10 * nrow(species)
}

city_metric_height = 30
traits_width = 50
phylo_tree_width = 125
title_height = 8

get_image_height = function(city_cluster_data_metrics) {
  get_presence_cell_height(city_cluster_data_metrics) + (3 * city_metric_height) + title_height
}

get_image_width = function(city_cluster_data_metrics) {
  get_presence_cell_width(city_cluster_data_metrics) + traits_width + phylo_tree_width
}

species_in_city_label = function(species) {
  paste(
    ifelse(species$seasonal == 'Resident', '', substring(species$seasonal, 0, 1)),
    ifelse(species$origin == 'Native', '', substring(species$origin, 0, 1)),
    ifelse(species$distance_to_northern_edge_km > 200, '', paste('NRL', round(species$distance_to_northern_edge_km), sep = ' ')),
    ifelse(species$distance_to_southern_edge_km > 200, '', paste('SRL', round(species$distance_to_northern_edge_km), sep = ' ')),
    sep = '\n'
  )
}

species_in_city_label(head(communities))

[1] "\nI\n\n" "\nI\n\n" "\nI\n\n" "P\n\n\n" "\n\n\n"  "\nI\n\n"

plot_city_cluster = function(city_cluster_data_metrics, title) {
  species_in_cluster = communities %>% 
    filter(city_id %in% city_cluster_data_metrics$city_id) %>% 
    dplyr::select(jetz_species_name, city_name, relative_abundance_proxy, seasonal, origin, distance_to_northern_edge_km, distance_to_southern_edge_km)
  
  species_in_cluster$label = species_in_city_label(species_in_cluster)
  tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
    
  gg_tree = ggtree(tree_cropped)
  
  gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) + 
          geom_tile(aes(fill=relative_abundance_proxy)) + 
          geom_text(aes(label=label), size=0.75) +
          scale_fill_gradientn(colours=c("#98FB98", "#FFFFE0", "yellow", "orange", "#FF4500", "red", "red"), values=c(0, 0.00000000001, 0.1, 0.25, 0.5, 0.75, 1), na.value = "transparent") +
          theme_minimal() + xlab(NULL) + ylab(NULL) + 
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
          labs(fill='Urban Proxy Abundance')
  
  species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
  
  gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
  
  gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
  
  gg_cities_gape_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = gape_width_fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("Gape") + ylim(0, 1)
  
  gg_cities_loco_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = handwing_index_fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("HWI") + ylim(0, 1)
  
  gg_title = ggplot() + labs(title = title, subtitle = paste(
    test_value('MNTD', city_cluster_data_metrics$mntd_normalised),
    test_value('Compound FDiv', city_cluster_data_metrics$fdiv_normalised),
    test_value('HWI FDiv', city_cluster_data_metrics$handwing_index_fdiv_normalised),
    test_value('Gape width FDiv', city_cluster_data_metrics$gape_width_fdiv_normalised),
    sep = '\n'
  )) + theme_minimal() + theme(plot.subtitle=element_text(size=8, hjust=0, color="#444444"))
  
  gg_presence_height = get_presence_cell_height(city_cluster_data_metrics)
  gg_presence_width = get_presence_cell_width(city_cluster_data_metrics)
  
  gg_presence %>% insert_top(gg_cities_loco_fd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_gape_fd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_mntd, height = (city_metric_height / gg_presence_height)) %>% insert_left(gg_tree, width = (phylo_tree_width / gg_presence_width)) %>% insert_right(gg_traits, width = (traits_width / gg_presence_width)) %>% insert_top(gg_title, height = (title_height / gg_presence_height))
}

REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive_using_abundance')

Nearctic

nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "San Jose"                 "Los Angeles"              "Concord"                  "Tijuana"                  "Bakersfield"             
 [6] "Fresno"                   "Sacramento"               "Mexicali"                 "Hermosillo"               "Las Vegas"               
[11] "Phoenix"                  "Tucson"                   "Durango"                  "Portland"                 "Chihuahua"               
[16] "Aguascalientes"           "Seattle"                  "Ciudad Juárez"            "San Luis Potosí"          "Mexico City"             
[21] "Saltillo"                 "Vancouver"                "Salt Lake City"           "Albuquerque"              "Monterrey"               
[26] "Nuevo Laredo"             "San Antonio"              "Denver"                   "Austin"                   "Houston"                 
[31] "Dallas"                   "Oklahoma City"            "Calgary"                  "New Orleans"              "Kansas City"             
[36] "Omaha"                    "St. Louis"                "Bradenton"                "Tampa"                    "Minneapolis [Saint Paul]"
[41] "Atlanta"                  "Orlando"                  "Louisville"               "Chicago"                  "Indianapolis"            
[46] "Milwaukee"               

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id)) 

Run 0 stress 0.1005678 
Run 1 stress 0.1000046 
... New best solution
... Procrustes: rmse 0.007231129  max resid 0.03470438 
Run 2 stress 0.1007176 
Run 3 stress 0.1017503 
Run 4 stress 0.1017503 
Run 5 stress 0.1217441 
Run 6 stress 0.1017503 
Run 7 stress 0.121519 
Run 8 stress 0.1210168 
Run 9 stress 0.1017503 
Run 10 stress 0.1000046 
... New best solution
... Procrustes: rmse 0.000008641127  max resid 0.00003146381 
... Similar to previous best
Run 11 stress 0.1217145 
Run 12 stress 0.1212382 
Run 13 stress 0.1323942 
Run 14 stress 0.102946 
Run 15 stress 0.1017503 
Run 16 stress 0.1017503 
Run 17 stress 0.1000046 
... Procrustes: rmse 0.000005711568  max resid 0.00001968977 
... Similar to previous best
Run 18 stress 0.1217145 
Run 19 stress 0.1005678 
Run 20 stress 0.1007176 
*** Best solution repeated 2 times

nearctic_cities_nmds

scree_plot(nearctic_cities_nmds)

nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 4)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(nearctic_cities)

nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = nearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Nearctic MNTD', nearctic_cities$mntd_normalised)

Warning: cannot compute exact p-value with ties

[1] "Nearctic MNTD mean 0.6 "

test_value('Nearctic beak gape FDiv', nearctic_cities$gape_width_fdiv_normalised)

Warning: cannot compute exact p-value with ties

[1] "Nearctic beak gape FDiv mean 0.59 "

test_value('Nearctic HWI FDiv', nearctic_cities$handwing_index_fdiv_normalised)

Warning: cannot compute exact p-value with ties

[1] "Nearctic HWI FDiv mean 0.31 **"

nrow(nearctic_cities)

[1] 46

Neartic Cluster 1`

nearactic_cluster1 = nearctic_cities %>% filter(cluster == 1)
plot_city_cluster(nearactic_cluster1, 'Neartic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'), width = get_image_width(nearactic_cluster1), height = get_image_height(nearactic_cluster1), units = "mm")

Neartic Cluster 2

nearactic_cluster2 = nearctic_cities %>% filter(cluster == 2)
plot_city_cluster(nearactic_cluster2, 'Neartic cluster 2')

Warning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'), width = get_image_width(nearactic_cluster2), height = get_image_height(nearactic_cluster2), units = "mm")

Neartic Cluster 3

nearactic_cluster3 = nearctic_cities %>% filter(cluster == 3)
plot_city_cluster(nearactic_cluster3, 'Neartic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'), width = get_image_width(nearactic_cluster3), height = get_image_height(nearactic_cluster3), units = "mm")

Neartic Cluster 4

nearactic_cluster4 = nearctic_cities %>% filter(cluster == 4)
plot_city_cluster(nearactic_cluster4, 'Neartic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'), width = get_image_width(nearactic_cluster4), height = get_image_height(nearactic_cluster4), units = "mm")

Neotropic

neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Culiacán"                  "Guadalajara"               "Morelia"                   "Acapulco"                  "Querétaro"                
 [6] "Cuernavaca"                "Puebla"                    "Oaxaca"                    "Xalapa"                    "Veracruz"                 
[11] "Tuxtla Gutiérrez"          "Villahermosa"              "Guatemala City"            "San Salvador"              "San Pedro Sula"           
[16] "Mérida"                    "Tegucigalpa"               "Managua"                   "San José"                  "Cancún"                   
[21] "Guayaquil"                 "Chiclayo"                  "Panama City"               "Trujillo"                  "Quito"                    
[26] "Havana"                    "Cali"                      "Lima"                      "Pereira"                   "Miami"                    
[31] "Medellín"                  "Ibagué"                    "Cartagena"                 "Kingston"                  "Bogota"                   
[36] "Barranquilla"              "Bucaramanga"               "Cúcuta"                    "Maracaibo"                 "Arequipa"                 
[41] "Barquisimeto"              "Santo Domingo"             "Maracay"                   "El Alto [La Paz]"          "Caracas"                  
[46] "Cochabamba"                "Viña del Mar [Valparaíso]" "Río Piedras [San Juan]"    "Barcelona"                 "Concepción"               
[51] "Santiago"                  "Mendoza"                   "Salta"                     "Cordoba"                   "Asuncion"                 
[56] "Buenos Aires"              "La Plata"                  "Ciudad del Este"           "Montevideo"                "Mar del Plata"            
[61] "Porto Alegre"              "São Paulo"                 "Santos"                    "Sao Jose dos Campos"      

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id)) 

Run 0 stress 0.134619 
Run 1 stress 0.1405831 
Run 2 stress 0.1402343 
Run 3 stress 0.1346226 
... Procrustes: rmse 0.001163646  max resid 0.007319248 
... Similar to previous best
Run 4 stress 0.1414222 
Run 5 stress 0.134433 
... New best solution
... Procrustes: rmse 0.00682457  max resid 0.04571633 
Run 6 stress 0.1344509 
... Procrustes: rmse 0.001186994  max resid 0.006732139 
... Similar to previous best
Run 7 stress 0.134433 
... New best solution
... Procrustes: rmse 0.00001887778  max resid 0.00006224176 
... Similar to previous best
Run 8 stress 0.1406946 
Run 9 stress 0.1346225 
... Procrustes: rmse 0.007127633  max resid 0.04569628 
Run 10 stress 0.1346406 
... Procrustes: rmse 0.007281613  max resid 0.04569614 
Run 11 stress 0.134433 
... Procrustes: rmse 0.00009355504  max resid 0.0002972054 
... Similar to previous best
Run 12 stress 0.1405822 
Run 13 stress 0.1405822 
Run 14 stress 0.1346226 
... Procrustes: rmse 0.007186606  max resid 0.04568917 
Run 15 stress 0.1346226 
... Procrustes: rmse 0.007196645  max resid 0.04566584 
Run 16 stress 0.134619 
... Procrustes: rmse 0.0068466  max resid 0.04566077 
Run 17 stress 0.134433 
... New best solution
... Procrustes: rmse 0.00001770911  max resid 0.00005326707 
... Similar to previous best
Run 18 stress 0.1346369 
... Procrustes: rmse 0.006925503  max resid 0.04573761 
Run 19 stress 0.141222 
Run 20 stress 0.1348046 
... Procrustes: rmse 0.006551197  max resid 0.04607941 
*** Best solution repeated 1 times

neotropic_cities_nmds

scree_plot(neotropic_cities_nmds)

neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(neotropic_cities)

neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = neotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Neotropic MNTD', neotropic_cities$mntd_normalised)

[1] "Neotropic MNTD mean 0.53 "

test_value('Neotropic beak gape FDiv', neotropic_cities$gape_width_fdiv_normalised)

[1] "Neotropic beak gape FDiv mean 0.44 "

test_value('Neotropic HWI FDiv', neotropic_cities$handwing_index_fdiv_normalised)

[1] "Neotropic HWI FDiv mean 0.53 "

nrow(neotropic_cities)

[1] 64

Neotropic Cluster 1

neotropic_cluster1 = neotropic_cities %>% filter(cluster == 1)
plot_city_cluster(neotropic_cluster1, 'Neotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'), width = get_image_width(neotropic_cluster1), height = get_image_height(neotropic_cluster1), units = "mm")

Neotropic Cluster 2

neotropic_cluster2 = neotropic_cities %>% filter(cluster == 2)
plot_city_cluster(neotropic_cluster2, 'Neotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'), width = get_image_width(neotropic_cluster2), height = get_image_height(neotropic_cluster2), units = "mm")

Neotropic Cluster 3

neotropic_cluster3 = neotropic_cities %>% filter(cluster == 3)
plot_city_cluster(neotropic_cluster3, 'Neotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'), width = get_image_width(neotropic_cluster3), height = get_image_height(neotropic_cluster3), units = "mm")

Neotropic Cluster 4

neotropic_cluster4 = neotropic_cities %>% filter(cluster == 4)
plot_city_cluster(neotropic_cluster4, 'Neotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'), width = get_image_width(neotropic_cluster4), height = get_image_height(neotropic_cluster4), units = "mm")

Neotropic Cluster 5

neotropic_cluster5 = neotropic_cities %>% filter(cluster == 5)
plot_city_cluster(neotropic_cluster5, 'Neotropic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'), width = get_image_width(neotropic_cluster5), height = get_image_height(neotropic_cluster5), units = "mm")

Palearctic

palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Lisbon"                "Porto"                 "Marrakesh"             "Seville"               "Dublin"                "Málaga"               
 [7] "Madrid"                "Glasgow"               "Bilbao"                "Liverpool"             "Bristol"               "Manchester"           
[13] "Birmingham"            "Leeds"                 "Newcastle upon Tyne"   "Sheffield"             "Nottingham"            "Valencia"             
[19] "London"                "Toulouse"              "Paris"                 "Barcelona"             "Rotterdam [The Hague]" "Brussels"             
[25] "Amsterdam"             "Lyon"                  "Marseille"             "Dusseldorf"            "Nice"                  "Frankfurt am Main"    
[31] "Zurich"                "Oslo"                  "Stuttgart"             "Hamburg"               "Genoa"                 "Nuremberg"            
[37] "Copenhagen"            "Munich"                "Berlin"                "Dresden"               "Rome"                  "Prague"               
[43] "Stockholm"             "Poznan"                "Vienna"                "Wroclaw"               "Zagreb"                "Gdansk"               
[49] "Budapest"              "Krakow"                "Warsaw"                "Helsinki"              "Riga"                  "Belgrade"             
[55] "Lviv"                  "Sofia"                 "Thessaloniki"          "Saint Petersburg"      "Minsk"                 "Athens"               
[61] "Kyiv"                  "Istanbul"              "Odesa"                 "Samsun"                "Luxor"                 "Tel Aviv"             
[67] "Jerusalem"             "Tbilisi"               "Yerevan"               "Kuwait City"           "Doha"                  "Abu Dhabi"            
[73] "Dubai"                 "Bishkek"              

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id)) 

Run 0 stress 0.04961857 
Run 1 stress 0.07367425 
Run 2 stress 0.05453968 
Run 3 stress 0.0712512 
Run 4 stress 0.05826357 
Run 5 stress 0.06338732 
Run 6 stress 0.06304725 
Run 7 stress 0.08618201 
Run 8 stress 0.07925636 
Run 9 stress 0.05647129 
Run 10 stress 0.09129792 
Run 11 stress 0.0497659 
... Procrustes: rmse 0.005703935  max resid 0.01420466 
Run 12 stress 0.05000939 
... Procrustes: rmse 0.06743721  max resid 0.2226243 
Run 13 stress 0.07603764 
Run 14 stress 0.0500095 
... Procrustes: rmse 0.06741103  max resid 0.2225553 
Run 15 stress 0.096573 
Run 16 stress 0.0556821 
Run 17 stress 0.05115583 
Run 18 stress 0.08625966 
Run 19 stress 0.06766821 
Run 20 stress 0.06829368 
Run 21 stress 0.07723847 
Run 22 stress 0.05236354 
Run 23 stress 0.05642252 
Run 24 stress 0.05235961 
Run 25 stress 0.06829346 
Run 26 stress 0.05353311 
Run 27 stress 0.06421851 
Run 28 stress 0.05706493 
Run 29 stress 0.0512458 
Run 30 stress 0.05479168 
*** Best solution was not repeated -- monoMDS stopping criteria:
    14: no. of iterations >= maxit
    14: stress ratio > sratmax
     2: scale factor of the gradient < sfgrmin

palearctic_cities_nmds

scree_plot(palearctic_cities_nmds)

Warning: Quick-TRANSfer stage steps exceeded maximum (= 3700)

palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 7)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(palearctic_cities)

palearctic_biomes = st_crop(resolve[resolve$REALM == 'Palearctic',c('REALM')], xmin = -30, ymin = 20, xmax = 80, ymax = 65)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = palearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Palearctic MNTD', palearctic_cities$mntd_normalised)

[1] "Palearctic MNTD mean 0.58 *"

test_value('Palearctic beak gape FDiv', palearctic_cities$gape_width_fdiv_normalised)

[1] "Palearctic beak gape FDiv mean 0.86 ***"

test_value('Palearctic HWI FDiv', palearctic_cities$handwing_index_fdiv_normalised)

[1] "Palearctic HWI FDiv mean 0.51 "

nrow(palearctic_cities)

[1] 74

Palearctic Cluster 1

palearctic_cluster1 = palearctic_cities %>% filter(cluster == 1)

plot_city_cluster(palearctic_cluster1, 'Palearctic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'), width = get_image_width(palearctic_cluster1), height = get_image_height(palearctic_cluster1), units = "mm")

Palearctic Cluster 2

palearctic_cluster2 = palearctic_cities %>% filter(cluster == 2)
plot_city_cluster(palearctic_cluster2, 'Palearctic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'), width = get_image_width(palearctic_cluster2), height = get_image_height(palearctic_cluster2), units = "mm")

Palearctic Cluster 3

palearctic_cluster3 = palearctic_cities %>% filter(cluster == 3)
plot_city_cluster(palearctic_cluster3, 'Palearctic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'), width = get_image_width(palearctic_cluster3), height = get_image_height(palearctic_cluster3), units = "mm")

Palearctic Cluster 4

palearctic_cluster4 = palearctic_cities %>% filter(cluster == 4)
plot_city_cluster(palearctic_cluster4, 'Palearctic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'), width = get_image_width(palearctic_cluster4), height = get_image_height(palearctic_cluster4), units = "mm")

Palearctic Cluster 5

palearctic_cluster5 = palearctic_cities %>% filter(cluster == 5)
nrow(palearctic_cluster5)

[1] 50

plot_city_cluster(palearctic_cluster5[1:25,], 'Palearctic cluster 5a')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'), width = get_image_width(palearctic_cluster5), height = get_image_height(palearctic_cluster5), units = "mm")

plot_city_cluster(palearctic_cluster5[25:50,], 'Palearctic cluster 5b')

Warning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5b.jpg'), width = get_image_width(palearctic_cluster5), height = get_image_height(palearctic_cluster5), units = "mm")

Palearctic Cluster 6

palearctic_cluster6 = palearctic_cities %>% filter(cluster == 6)
plot_city_cluster(palearctic_cluster6, 'Palearctic cluster 6')

Warning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'), width = get_image_width(palearctic_cluster6), height = get_image_height(palearctic_cluster6), units = "mm")

Palearctic Cluster 7

palearctic_cluster7 = palearctic_cities %>% filter(cluster == 7)
plot_city_cluster(palearctic_cluster7, 'Palearctic cluster 7')

Warning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster7.jpg'), width = get_image_width(palearctic_cluster7), height = get_image_height(palearctic_cluster7), units = "mm")

Afrotropic

afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
[1] "Cape Town"    "Johannesburg" "Pretoria"     "Kigali"       "Kampala"      "Arusha"       "Nairobi"      "Addis Ababa"  "Antananarivo"

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id)) 

Run 0 stress 0.00009014786 
Run 1 stress 0.0007443902 
Run 2 stress 0.000377487 
... Procrustes: rmse 0.002074146  max resid 0.003340817 
... Similar to previous best
Run 3 stress 0.00008179516 
... New best solution
... Procrustes: rmse 0.0001660244  max resid 0.000361027 
... Similar to previous best
Run 4 stress 0.000656661 
Run 5 stress 0.3209238 
Run 6 stress 0.00009538562 
... Procrustes: rmse 0.00003771484  max resid 0.00007382498 
... Similar to previous best
Run 7 stress 0.00009646529 
... Procrustes: rmse 0.0002653617  max resid 0.0004782319 
... Similar to previous best
Run 8 stress 0.00009486991 
... Procrustes: rmse 0.0002901629  max resid 0.0004478019 
... Similar to previous best
Run 9 stress 0.00009857688 
... Procrustes: rmse 0.0004629411  max resid 0.0009418255 
... Similar to previous best
Run 10 stress 0.00008547616 
... Procrustes: rmse 0.00001297639  max resid 0.00002426197 
... Similar to previous best
Run 11 stress 0.00009878898 
... Procrustes: rmse 0.0004635066  max resid 0.0009419129 
... Similar to previous best
Run 12 stress 0.0005427305 
... Procrustes: rmse 0.002940686  max resid 0.004342951 
... Similar to previous best
Run 13 stress 0.001227374 
Run 14 stress 0.00009793901 
... Procrustes: rmse 0.0001331878  max resid 0.0002796824 
... Similar to previous best
Run 15 stress 0.00009902405 
... Procrustes: rmse 0.0001283685  max resid 0.0002454039 
... Similar to previous best
Run 16 stress 0.00009120045 
... Procrustes: rmse 0.00002743712  max resid 0.00005350613 
... Similar to previous best
Run 17 stress 0.00009094513 
... Procrustes: rmse 0.0002505202  max resid 0.00046142 
... Similar to previous best
Run 18 stress 0.00009137309 
... Procrustes: rmse 0.000123082  max resid 0.0002401071 
... Similar to previous best
Run 19 stress 0.00009536048 
... Procrustes: rmse 0.00003982921  max resid 0.00006129269 
... Similar to previous best
Run 20 stress 0.00006560657 
... New best solution
... Procrustes: rmse 0.0001813766  max resid 0.0003094943 
... Similar to previous best
*** Best solution repeated 1 times

Warning: stress is (nearly) zero: you may have insufficient data

afrotropic_cities_nmds

scree_plot(afrotropic_cities_nmds)

afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 2)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(afrotropic_cities)

afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Afrotropic MNTD', afrotropic_cities$mntd_normalised)

[1] "Afrotropic MNTD mean 0.14 *"

test_value('Afrotropic beak gape FDiv', afrotropic_cities$gape_width_fdiv_normalised)

[1] "Afrotropic beak gape FDiv mean 0.41 "

test_value('Afrotropic HWI FDiv', afrotropic_cities$handwing_index_fdiv_normalised)

[1] "Afrotropic HWI FDiv mean 0.61 "

nrow(afrotropic_cities)

[1] 9

Afrotropic Cluster 1

afrotropic_cluster1 = afrotropic_cities %>% filter(cluster == 1)
plot_city_cluster(afrotropic_cluster1, 'Afrotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'), width = get_image_width(afrotropic_cluster1), height = get_image_height(afrotropic_cluster1), units = "mm")

Afrotropic Cluster 2

afrotropic_cluster2 = afrotropic_cities %>% filter(cluster == 2)
plot_city_cluster(afrotropic_cluster2, 'Afrotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'), width = get_image_width(afrotropic_cluster2), height = get_image_height(afrotropic_cluster2), units = "mm")

Indomalayan

indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
  [1] "Srinagar"            "Jamnagar"            "Jammu"               "Rajkot"              "Bikaner"             "Jodhpur"             "Jalandhar"          
  [8] "Ahmedabad"           "Bhavnagar"           "Ludhiana"            "Anand"               "Udaipur"             "Surat"               "Vadodara"           
 [15] "Ajmer"               "Chandigarh"          "Vasai-Virar"         "Mumbai"              "Jaipur"              "Delhi [New Delhi]"   "Nashik"             
 [22] "Dehradun"            "Kota"                "Pune"                "Haridwar"            "Dhule"               "Ujjain"              "Indore"             
 [29] "Ahmadnagar"          "Kolhapur"            "Jalgaon"             "Agra"                "Aurangabad"          "Sangli"              "Belagavi"           
 [36] "Gwalior"             "Budaun"              "Bareilly"            "Dharwad"             "Bhopal"              "Bhind"               "Mangaluru"          
 [43] "Solapur"             "Vijayapura"          "Akola"               "Latur"               "Kannur"              "Davanagere"          "Thalassery"         
 [50] "Amravati"            "Kalaburagi"          "Kozhikode"           "Guruvayur"           "Malappuram"          "Lucknow"             "Thrissur"           
 [57] "Mysuru"              "Kochi"               "Alappuzha"           "Nagpur"              "Kollam"              "Jabalpur"            "Ettumanoor"         
 [64] "Hyderabad"           "Coimbatore"          "Bengaluru"           "Thiruvananthapuram"  "Tiruppur"            "Faizabad"            "Erode"              
 [71] "Prayagraj"           "Pratapgarh"          "Salem"               "Dindigul"            "Madurai"             "Tiruchirappalli"     "Durg"               
 [78] "Vellore"             "Tirupati"            "Raipur"              "Bilaspur"            "Vijayawada"          "Puducherry"          "Chennai"            
 [85] "Kathmandu"           "Colombo"             "Rajamahendravaram"   "Patna"               "Kandy"               "Bihar Sharif"        "Visakhapatnam"      
 [92] "Ranchi"              "Brahmapur"           "Jamshedpur"          "Darjeeling"          "Siliguri"            "Cuttack"             "Bhubaneshwar"       
 [99] "Jalpaiguri"          "Berhampore"          "Kolkata"             "Krishnanagar"        "Guwahati [Dispur]"   "Agartala"            "Silchar"            
[106] "Dimapur"             "Bangkok"             "George Town"         "Kuala Lumpur"        "Phnom Penh"          "Singapore"           "Hong Kong"          
[113] "Sha Tin"             "Hsinchu"             "Taichung"            "New Taipei [Taipei]" "Tainan"              "Denpasar"            "Kaohsiung"          
[120] "Kota Kinabalu"      

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id)) 

Run 0 stress 0.1190668 
Run 1 stress 0.1224401 
Run 2 stress 0.1394961 
Run 3 stress 0.1161598 
... New best solution
... Procrustes: rmse 0.02798297  max resid 0.2346031 
Run 4 stress 0.1540252 
Run 5 stress 0.1175538 
Run 6 stress 0.1167657 
Run 7 stress 0.1384412 
Run 8 stress 0.1241275 
Run 9 stress 0.1163498 
... Procrustes: rmse 0.025336  max resid 0.2481842 
Run 10 stress 0.1153501 
... New best solution
... Procrustes: rmse 0.008287738  max resid 0.08217213 
Run 11 stress 0.1199235 
Run 12 stress 0.1172201 
Run 13 stress 0.151756 
Run 14 stress 0.1246307 
Run 15 stress 0.134023 
Run 16 stress 0.1580175 
Run 17 stress 0.1229198 
Run 18 stress 0.1502633 
Run 19 stress 0.1189366 
Run 20 stress 0.117072 
Run 21 stress 0.1458634 
Run 22 stress 0.137476 
Run 23 stress 0.1186188 
Run 24 stress 0.1637221 
Run 25 stress 0.1191928 
Run 26 stress 0.1164559 
Run 27 stress 0.142366 
Run 28 stress 0.1379337 
Run 29 stress 0.1191853 
Run 30 stress 0.1266232 
*** Best solution was not repeated -- monoMDS stopping criteria:
    29: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

indomalayan_cities_nmds

scree_plot(indomalayan_cities_nmds)

indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(indomalayan_cities)

indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
 
ggplot() + 
  geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) + 
  geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Indomalayan MNTD', indomalayan_cities$mntd_normalised)

[1] "Indomalayan MNTD mean 0.44 ***"

test_value('Indomalayan beak gape FDiv', indomalayan_cities$gape_width_fdiv_normalised)

[1] "Indomalayan beak gape FDiv mean 0.52 "

test_value('Indomalayan HWI FDiv', indomalayan_cities$handwing_index_fdiv_normalised)

[1] "Indomalayan HWI FDiv mean 0.82 ***"

nrow(indomalayan_cities)

[1] 120

Indomalayan Cluster 1

indomalayan_cluster1 = indomalayan_cities %>% filter(cluster == 1)
plot_city_cluster(indomalayan_cluster1, 'Indomalayan cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'), width = get_image_width(indomalayan_cluster1), height = get_image_height(indomalayan_cluster1), units = "mm")

Indomalayan Cluster 2

indomalayan_cluster2 = indomalayan_cities %>% filter(cluster == 2)
plot_city_cluster(indomalayan_cluster2, 'Indomalayan cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'), width = get_image_width(indomalayan_cluster2), height = get_image_height(indomalayan_cluster2), units = "mm")

Indomalayan Cluster 3

indomalayan_cluster3 = indomalayan_cities %>% filter(cluster == 3)
plot_city_cluster(indomalayan_cluster3, 'Indomalayan cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'), width = get_image_width(indomalayan_cluster3), height = get_image_height(indomalayan_cluster3), units = "mm")

Indomalayan Cluster 4

indomalayan_cluster4 = indomalayan_cities %>% filter(cluster == 4)
plot_city_cluster(indomalayan_cluster4, 'Indomalayan cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'), width = get_image_width(indomalayan_cluster4), height = get_image_height(indomalayan_cluster4), units = "mm")

Indomalayan Cluster 5

indomalayan_cluster5 = indomalayan_cities %>% filter(cluster == 5)
plot_city_cluster(indomalayan_cluster5, 'Indomalayan cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'), width = get_image_width(indomalayan_cluster5), height = get_image_height(indomalayan_cluster5), units = "mm")

Australasia

australasia_cities_community_data = community_data_metrics %>% filter(core_realm == 'Australasia')
australasia_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
[1] "Perth"     "Adelaide"  "Melbourne" "Sydney"    "Brisbane"  "Auckland" 

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

australasia_cities_nmds = community_nmds(communities %>% filter(city_id %in% australasia_cities_community_data$city_id)) 

Run 0 stress 0 
Run 1 stress 0 
... Procrustes: rmse 0.1120339  max resid 0.1622729 
Run 2 stress 0.1572548 
Run 3 stress 0 
... Procrustes: rmse 0.07245479  max resid 0.1129183 
Run 4 stress 0 
... Procrustes: rmse 0.04603582  max resid 0.07090683 
Run 5 stress 0 
... Procrustes: rmse 0.2182675  max resid 0.295989 
Run 6 stress 0 
... Procrustes: rmse 0.133075  max resid 0.2072423 
Run 7 stress 0.0000995744 
... Procrustes: rmse 0.1149167  max resid 0.1742941 
Run 8 stress 0 
... Procrustes: rmse 0.1694089  max resid 0.2635754 
Run 9 stress 0 
... Procrustes: rmse 0.1063225  max resid 0.1577545 
Run 10 stress 0 
... Procrustes: rmse 0.1313955  max resid 0.2135724 
Run 11 stress 0 
... Procrustes: rmse 0.1251736  max resid 0.1780552 
Run 12 stress 0 
... Procrustes: rmse 0.05487419  max resid 0.08977988 
Run 13 stress 0.1572548 
Run 14 stress 0.00006459844 
... Procrustes: rmse 0.1484605  max resid 0.2073383 
Run 15 stress 0.00008456015 
... Procrustes: rmse 0.1937182  max resid 0.3430883 
Run 16 stress 0 
... Procrustes: rmse 0.09381547  max resid 0.1446598 
Run 17 stress 0.1572548 
Run 18 stress 0.00001782659 
... Procrustes: rmse 0.132877  max resid 0.2076538 
Run 19 stress 0 
... Procrustes: rmse 0.115211  max resid 0.1946323 
Run 20 stress 0 
... Procrustes: rmse 0.1542908  max resid 0.2134235 
Run 21 stress 0.1572548 
Run 22 stress 0 
... Procrustes: rmse 0.1239607  max resid 0.190374 
Run 23 stress 0.00007261303 
... Procrustes: rmse 0.123741  max resid 0.1717834 
Run 24 stress 0.00006062953 
... Procrustes: rmse 0.1354649  max resid 0.2269844 
Run 25 stress 0.00008963805 
... Procrustes: rmse 0.1180163  max resid 0.1705588 
Run 26 stress 0.000009458018 
... Procrustes: rmse 0.1301844  max resid 0.1953992 
Run 27 stress 0.2349416 
Run 28 stress 0.0000881637 
... Procrustes: rmse 0.07152338  max resid 0.1085942 
Run 29 stress 0 
... Procrustes: rmse 0.09917667  max resid 0.1266893 
Run 30 stress 0 
... Procrustes: rmse 0.08881426  max resid 0.1290499 
*** Best solution was not repeated -- monoMDS stopping criteria:
    25: stress < smin
     1: stress ratio > sratmax
     4: scale factor of the gradient < sfgrmin

Warning: stress is (nearly) zero: you may have insufficient data

australasia_cities_nmds

scree_plot(australasia_cities_nmds)

australasia_cities = cluster_cities(city_nmds = australasia_cities_nmds, cities_community_data = australasia_cities_community_data, centers = 3)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(australasia_cities)

australasia_biomes = st_crop(resolve[resolve$REALM == 'Australasia',c('REALM')], xmin = 80, ymin = -70, xmax = 180, ymax = 20)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = australasia_biomes, aes(geometry = geometry)) + 
  geom_sf(data = australasia_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Australasia MNTD', australasia_cities$mntd_normalised)

[1] "Australasia MNTD mean 0.48 "

test_value('Australasia beak gape FDiv', australasia_cities$gape_width_fdiv_normalised)

[1] "Australasia beak gape FDiv mean 0.5 "

test_value('Australasia HWI FDiv', australasia_cities$handwing_index_fdiv_normalised)

[1] "Australasia HWI FDiv mean 0.82 "

nrow(australasia_cities)

[1] 6

Australasia Cluster 1

australasia_cluster1 = australasia_cities %>% filter(cluster == 1)
plot_city_cluster(australasia_cluster1, 'Australasia cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_cluster1.jpg'), width = get_image_width(australasia_cluster1), height = get_image_height(australasia_cluster1), units = "mm")

australasia_cluster2 = australasia_cities %>% filter(cluster == 2)
plot_city_cluster(australasia_cluster2, 'Australasia cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_cluster2.jpg'), width = get_image_width(australasia_cluster2), height = get_image_height(australasia_cluster2), units = "mm")

australasia_cluster3 = australasia_cities %>% filter(cluster == 3)
plot_city_cluster(australasia_cluster3, 'Australasia cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_cluster3.jpg'), width = get_image_width(australasia_cluster3), height = get_image_height(australasia_cluster3), units = "mm")